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We analyze the quantum dynamics of two electromagnetic oscillators coupled in series to a voltage 
biased Josephson junction. When the applied voltage leads to a Josephson frequency across the 
junction which matches the sum of the two mode frequencies, tunneling Cooper pairs excite photons 
in both modes simultaneously leading to far-from-equilibrium states. These states display highly 
non-classical features including strong anti-bunching, violation of Cauchy-Schwartz inequalitites, 
and number squeezing. The regimes of low and high photon occupancies allow for analytical results 
which are supported by a full numerical treatment. The impact of asymmetries between the two 
modes is explored, revealing a pronounced enhancement of number squeezing when the modes are 
damped at different rates. 
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I. INTRODUCTION 

It has long been known that the current flowing 
through a voltage-biased mesoscopic conductor can pro¬ 
vide an extremely sensitive probe of its electromagnetic 
environment— ~—. The current-voltage characteristics of 
a tunnel junction placed in series with a transmission 
line resonator is a particularly well-studied case— 
The transmission line resonator contains a series of well- 
defined harmonic inodes whose presence opens up inelas¬ 
tic current channels leading to characteristic features in 
the dc current flowing through the junction—. The ad¬ 
vent of high-Q superconducting resonators whose quan¬ 
tum state can be measured with great precision— together 
with the development of hybrid devices which couple non- 
metallic conductors to resonators—£, has led to a renewed 
interest in the interaction between tunneling electrons or 
Cooper pairs and harmonic modes. Whilst earlier ex¬ 
periments— ^ on mesoscopic conductors coupled to elec¬ 
tromagnetic resonators focussed on how the harmonic 
modes affect the current in a regime where the modes 
themselves are close to thermal equilibrium, more recent 
experimental— - — and theoretical work— — has begun 
to investigate how the current influences the resonator 
state and to explore the dynamics of systems where the 
resonator is far from thermal equilibrium. 

For a Josephson junction which is biased with a sub¬ 
gap voltage, V , the relationship between the dc current 
and the energy pumped into the electromagnetic environ¬ 
ment is particularly simple as all of the energy associated 
with a tunneling Cooper pair must be absorbed by the 
environment jii. When the Josephson junction is placed 
in series with a transmission line resonator a dc current is 
expected when the ac Josephson frequency toj = 2 eV/h 
matches one or more of the mode frequencies in the trans¬ 
mission line. Experiments using low-Q resonators— 
have demonstrated that when the individual harmonic 
modes remain close to thermal equilibrium, they lead to 
well-defined peaks in the dc current whose heights and 
widths can be calculated using perturbation theory. In 
contrast, a high-Q resonator can be excited to far-from- 
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FIG. 1: Effective circuit model of the system. It consists 
of a Josephson junction (J J) in series with two LC oscil¬ 
lators, across which a voltage V is applied. The two LC 
oscillators are assumed to have different angular frequencies 
uj a = (L a C a )~ 1/2 ±w b = (L b C b )~ 1/2 . 

equilibrium states containing many photons— which are 
predicted to display intriguing non-classical features such 
as number squeezing——. This new field of Josephson 
photonics combines typical processes known from quan¬ 
tum optical set-ups with those known from charge trans¬ 
fer physics in highly versatile devices. 

In this article we consider a voltage-biased supercon¬ 
ducting junction whose ac Josephson frequency is tuned 
to excite two electromagnetic modes simultaneously (see 
Fig. ED. Signatures of such processes have been observed 
in the dc current flowing through Josephson junctions 
coupled to low-Q resonators and can also be understood 
within a perturbative approach as the modes remain close 
to thermal equilibrium. While we address this domain as 
well, our main focus here lies in the regime where the 
power transferred to the resonator modes is sufficient to 
drive them into far-from-equilibrium states while still dis¬ 
playing strong quantum properties. Note that the system 
we consider here differs from those used in recent exper¬ 
iments to produce photon pairs 25 i 26 in that the energy 
comes from a dc voltage. 

Starting from a simple model Hamiltonian which de¬ 
scribes the effect of the Cooper pairs on the oscillators 
through a highly non-linear ac drive at the Josephson fre¬ 
quency, we use a rotating wave approximation to derive 
an effective time-independent Hamiltonian which we use 
to analyse the quantum dynamics of the oscillators. Al¬ 
though the full behavior of the system can only be uncov- 
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ered by numerical solutions of the quantum master equa¬ 
tion, we find that approximate analytical descriptions are 
available in the two regimes of low and high photon occu¬ 
pancy. In the former one a perturbative treatment in the 
Josephson energy applies while in the latter explicit re¬ 
sults are obtained by linearizing about the classical fixed 
points which provide a faithful description of the quan¬ 
tum dynamics when the zero-point fluctuations of the 
oscillators are small. 

The excitation of the two oscillators shows a clear 
threshold as a function of the Cooper pair pumping rate. 
Earlier work, which investigated the quantum dynamics 
of a single mode— — driven by a voltage-biased Joseph¬ 
son junction, showed that non-classical features in the 
state of the oscillator such as number squeezing (sub- 
Poissonian photon statistics) occur very generally. For 
the two-mode system, we also find significant number¬ 
squeezing occurs in the states of the individual oscilla¬ 
tors, especially in the above-threshold regime where the 
oscillators are strongly excited. Interestingly, when the 
damping rates of the oscillators are very unequal, the less- 
damped oscillator displays much stronger strong number 
squeezing than is ever found for a single-oscillator sys¬ 
tem. Provided that the quantum zero-point fluctuations 
are not too small, the number squeezing is strong enough 
to lead to negative regions in the Wigner function. 

This work is organised as follows. We introduce our 
theoretical model in Sec. eh and analyse its low photon 
limit in Sec. EH and its semi-classical dynamics in Sec. 
El Sections El and IVII explore the quantum dynamics 
of the system in the below and above threshold regimes, 
respectively. Finally, Sec. IVIII contains a discussion and 
the conclusions. The Appendix contains further details 
on some of the calculations described in the main text. 


Here we analyze the case where the system is oper¬ 
ated close to the resonance that occurs when the voltage 
energy lost by a single Cooper pair traversing the circuit 
matches the energy required to simultaneously create one 
photon in each of the LC oscillators, uj = 2 eV/h = 
uj a + tub- We assume that the modes are not degener¬ 
ate so that u) a ^ uib- This means that the resonance at 
uj = oj a + LOb does not compete with processes in which 
two photons are absorbed by just one of the modes. 

We examine the behavior of the system as a function of 
the Josephson energy which describes the strength of the 
Cooper pair tunneling. The value of Ej can be thought 
of like a pumping rate for the oscillators: as it is increased 
the oscillators will be more strongly driven, become more 
strongly excited and behave more non-linearly. In prac¬ 
tice Ej can be varied in an effective single-junction by 
forming two junctions in parallel and applying a tunable 
flux in the SQUID loop that they form 20 i 29 . 

The strengths of the quantum fluctuations parame- 
terised by A a , A b , also play a very interesting role in 
determining the dynamics of the system and we will ex¬ 
amine how the behavior is modified when they are varied. 
For systems where a Josephson junction is embedded in a 
superconducting resonator designed to have a very high- 
Q the quantum fluctuations will typically be very small 
A 0 (b) < 1 . However, significantly stronger quantum 
fluctuations have very recently been engineered in low-Q 
resonators coupled to tunnel junctions— and it may be 
possible to combine stronger quantum fluctuations with 
higher Q values in the future. 


A. Rotating wave approximation 


II. MODEL SYSTEM 

We consider a system consisting of a Josephson junc¬ 
tion in series with two LC oscillators, A and B with 
angular frequencies ui a and co b across which a voltage V 
is applied (see Fig. [[]). The two oscillators could both 
be modes of a single superconducting resonator in which 
a Josephson junction is embedded between the ground 
plane and center conductor 13 i 20 i 27 i 28 (See Ref. [20j for a 
detailed derivation of the Hamiltonian for this case), but 
the system could also be realized using modes of two dif¬ 
ferent electrical resonators—. The effective Hamiltonian 
of the system takes the form 

H = HuiadJa + fkv b b'b (1) 

-Ej cos [uijt + A a (a + a 1 ') + A b (b + 6 1 ')] , 

where Ej is the Josephson energy of the junction, a and 
b are the lowering operators of the oscillators with fre¬ 
quencies to a and Wb respectively, and uj = 2 eV/h. The 
parameters A a ( b ) quantify the strength of the zero-point 
fluctuations of the oscillators, A a (6) = (2e 2 Z a{b) /h) 1 / 2 
where Z a{b) = yjL a ( b )/C a ( b ) is the impedance. 


The explicit time-dependence in the Hamiltonian com¬ 
plicates the analysis of the corresponding dynamics sig¬ 
nificantly. However, close to the resonance we are in¬ 
terested in, ujj ~ uj a + ui b , only some of the terms will 
play an important role and these can be picked out by a 
rotating wave approximation (RWA). 

We proceed following the approach in Refs. l20l - [22l We 
move to a rotating frame, applying a unitary transforma¬ 
tion of the form U{t) = e *w a a t at e iw i) 6 t 6t w } iere we define 
uj a +oj b — ivj, and make a RWA in which we neglect all of 
the rapidly oscillating terms in the rotating frame. The 
resulting effective Hamiltonian takes the form, 


H rwa = hS^aia + hS^tfb 


( 2 ) 


+ 


Ej J 1 (2A a VaJa)J 1 (2A b Vtfb) , t t s 

— {ab +ab): ’ 


where the colons imply normal ordering of the op¬ 
erators, (jW = — (Jji for i = a, b and Ej = 

Uje _ ( A “ +A| >)/ 2 . For sufficie ntly low photon numbers 
(such that 2A Q i/(ata), 2A b \/(6f&) -C 1) we can expand 
the Bessel functions in Eq. m to lowest order. In this 
limit the system reduces to a non-degenerate parametric 
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amplifier— 

i4w A = M (a) « t a+M (b) b t b+ ^ jA ^ Ab ( a tftt + o6 ) . ( 3 ) 


B. Quantum master equation 

The two oscillators are assumed to be weakly damped 
at rates 7 a and 75 which in general will not be the same. 
We therefore assume that the quantum master equation 
of the system takes the standard quantum optical form 
in the T = 0 limit— 

^ = -j[Hrwa, p\ + ^ (2 apa* - a ] ap - pa}a) 

+Y~ (2 bpb^ — bp — pb^b) , (4) 

where we adopt dimensionless units of time r = ty/j a 7b, 

r = yfuhb and Hrwa = H RWA /(h y /'y a 7b)- 

In an actual experimental realization of the JJ- 
oscillators system in Fig. [T] the damping of the oscil¬ 
lators (due to photon decay from the resonators) is not 
the only source of dissipation. Indeed, the existence and 
impact of local voltage fluctuations at the JJ can be seen 
in the broadening of the spectrum of emitted microwave 
radiation 11 1 21 . The existence of such fluctuations neces¬ 
sitates including explicitly an extra degree of freedom 
for the number of Cooper pairs N transported across 
the junction in the model. In the effective Hamilto¬ 
nian, Eq. ©, the (aW + ab ) term then gets replaced by 
(e ZT) odbt + e~ lv ab ), where = J2n |.W)(iV±l|. Local 
voltage fluctuations are included by an additional dissi- 
pator in © which in the simplest version takes the form 
C\N p] = r j(2Np N— N 2 p— pN 2 ) with rj = 7j/y/la7b- 
Ref. [2ll describes, how to treat the corresponding quan¬ 
tum master equation in the extended JJ-resonator space. 

However, it turns out that only certain observables 
sensitively depend on the strength of these fluctuations, 
characterized by 7 j, for example the spectral broadening. 
For other observables, such as the photon occupation and 
photonic correlation functions that are of relevance for 
this work, the impact of local voltage fluctuations is ba¬ 
sically negligible since experimentally one typically has 
7 j -C 7 a ,b (see for example Ref. |ll|). Then, formally, 
the Hamiltonian m is regained by putting 7 j = 0 so 
that the phase operators e ±lv simply appear as phase 
factors which can be removed via the gauge transforma¬ 
tion e lT1 / 2 a\e lT1 / 2 b^ —> a^, 6 b Note that this reflects a 
phase invariance of the RWA Hamiltonian ©. 


C. Relevant observables 

The basic structure of the RWA Hamiltonian [Eq. ©] 
in which photons are always created (or destroyed) jointly 
in the two oscillators and the linear damping that we as¬ 
sumed in formulating the master equation lead to a sim¬ 
ple connection between the occupation numbers of the 


two modes n a ^) = {a'a(b'b)) and the average dc current, 
Idc, flowing through the junction that can be obtained 
from an energy balance argument without the need to 
work with a current operator. Since each Cooper-pair 
that contributes to the dc current must create exactly 
one additional photon in each of the oscillators, the re¬ 
quirement that the energy gain and loss rates balance 
tells us that 


= 7a n a = 7 bn b , (5) 

2 e 

where in this case we have returned to dimension-full 
units. 

The quantum nature of the photonic states in the os¬ 
cillators is captured by photon correlation functions such 
as 




Jatfb) 


a(b) 


n a n b 


and the Fano factors 


( 6 ) 


F a {b) — 


([a + a.( 6 t 6 )] 2 > - n 2 a(b) 


1a(b) 


( 7 ) 


Whilst these two types of correlation functions are closely 
related to each other, they are nevertheless useful to char¬ 
acterize the photonic states in opposite regimes of pa¬ 
rameter space. In the regime of weak driving and low 
photon occupation deviations from the case of a driven 
harmonic oscillator are best seen in the g^ functions. 
Namely, with increasing driving amplitude Ej , the pho¬ 
ton distributions for the number states in the cavities 
evolve from Poissonian distributions with almost empty 
cavities towards distributions peaked around finite mean 
occupations n a ,n b . In this case the g^(O) functions 
© sensitively indicate deviations from the linear regime 

9aa{bb) (°) = 1 with 9ab (°) ^ 0 capturing growing cavity- 
cavity correlations. In the opposite regime of strong driv¬ 
ing, nonlinearities may substantially influence the widths 
of the peaks for photon occupations (energy fluctuations) 
as properly measured in the Fano-factors ©■ 

In the following, we will first focus on the regime of 
low photon occupancy, where charge transfer through the 
Josephson junction occurs sequentially (Coulomb block¬ 
ade regime) and analytical results can be obtained via a 
perturbative treatment in the drive amplitude Ej. In 
the opposite domain of large photons numbers in the 
cavities, a semi-classical type of approximation applies, 
though, with substantial quantum fluctuations. In both 
domains, the magnitude of the parameters A a ,A b , i.e. 
the strength of the ground state fluctuations of the re¬ 
spective cavities, plays a decisive role: they rule the im¬ 
pact of nonlinearities in the former regime and control the 
impact of quantum effects in the latter one. Correspond¬ 
ing findings will be supported by numerical solutions of 
the stationary states according to m and ©■ 
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III. FEW-PHOTON LIMIT 


The physics of the system described by the Hamilto¬ 
nian © and the master equation © is at its simplest 
when it is driven so weakly that excitations in the res¬ 
onators will relax to equilibrium well before a new excita¬ 
tion occurs. In that regime, very few photons, n a /b <C 1, 
reside in the resonators on average. Transport across the 
junction in turn is in the (dynamical) Coulomb-blockade 
regime, where subsequent Cooper-pair tunneling events 
occur uncorrelated with some tunneling rate. While the 
charge flows uncorrelated, the photons exhibit correla¬ 
tions already at the weakest driving. 

Now, for the present set-up one derives from the full 
RWA Hamiltonian 


n a 


iEj 

2 E c jr 


{: (ab-a^tf) 


Ji( 2 A a 


' a'a 


Ji( 2A aV / frtfr) 

A bVb^b 

( 8 ) 


with Ej = (A v / 7 a 7 b/A a A h )e( A “ +A i >)/ 2 and where rib fol¬ 
lows by replacing r —> 1 jr. In the lowest order in the 
driving strength this reduces to 


(o) = I (E ±\ 2 _i±i!_ 

4 \Ej J r 2 (<5(“) +<5 ( fc )) 2 + (l + r 2 ) 2 /4’ 


(9) 


with the superscript indicating the leading order in E 2 
and with again following from r —> 1/r. 




FIG. 2: (Color online) Autocorrelations 3 ^( 6 i) )(0) (left) and 

cross-correlations 3 ^( 0 ) (right) of the two modes vary with 
the strength of zero-point fluctuations A a (b) in the two oscil¬ 
lators. For weak driving, Ej = 0.2 Ej, the autocorrelations 
(symbols) are given by (11211 (lines) when A a , or simultane¬ 
ously A a and Af, are tuned. The reduced cross-correlations 
S , a 6 > ( 0 ) — l/( 2 n) (lines) obey the general relation (Hill with the 
mean of the autocorrelations [g^S (0) + 3 m^( 0)]/2 depicted as 
symbols for the case of symmetric damping r = 1 . 

For the correlations we focus on the symmetric case 
7 a = 7 b at resonance so that n a = rib = n. Then, based 
on the full Hamiltonian © one can show the general 
relation 


(a^ab^b) 



( 10 ) 


which implies 


9T(0) = l + - 2 [ 9 m + 9i^)_ 

with n as given in ©. Now, working to order Ej, one 
finds 


1 r 


( 11 ) 


( 12 ) 

Two types of correlations are encoded in the above 
3 ^( 0 ) functions. The most obvious ones stem from the 
common excitation process of photons in the two res¬ 
onators. They are therefore already present in the para¬ 
metric amplifier limit of the Hamiltonian © and well 
understood for that case, see e.g. Ref. HH- A convenient 
tool to characterise them is the noise reduction factor— 
NRF = [((ala — b^b) 2 ) — (n a — rib) 2 }/{n a + rib ) which in 
the symmetric situation y a = yb takes the form 


NRF = - 
2 


<?i 2 J( 0 ) + g%\o) -2ffS ) (0) 


+ 1 . 


(13) 


However, the perfect correlation of the excitation pro¬ 
cess leads to perfectly correlated occupations in the os¬ 
cillators with a noise reduction factor NRF = 0 only for 
the undamped case y Q = 7 b = 0. For any finite photon 
lifetimes in the cavities, the decay out of the two cavi¬ 
ties occurs uncorrelated which according to HID always 
implies in the stationary state and for the symmetric sit¬ 
uation NRF = 1/2. 

Further correlations in the light field are caused by the 
back-action of the resonator occupations on the photon 
creation processes. Generally speaking, the existence of 
photonic excitations in the resonators can either increase 
the probability of further excitations, similar to a stim¬ 
ulated emission effect, or it can hinder further excita¬ 
tions. Formally, these effects are encoded in the tran¬ 
sition matrix elements of the RWA-Hamiltonian © be¬ 
tween neighboring oscillator states, where the nonlinear¬ 
ities of the Bessel functions enter. If charge quantization 
of the Cooper-pair current is significant, the parameters 
A a /b become large, so that the nonlinearities already ap¬ 
pear at the few photon level. For the case of a single 
resonator, it was shown in Ref. [2l| that A 2 = 2 can com¬ 
pletely suppress transitions to higher occupations and re¬ 
duces the resonator effectively to a two-level system, thus 
operating as a perfect single photon source. The behav¬ 
ior of the correlation functions in the two-mode case is 
shown in Fig. [2] While a non-zero g'if,! (0) requires oscil¬ 
lator a to be populated up to the second excited state by 
two successive photons, this need not be the case for os¬ 
cillator b as it can relax before the second photon arrives. 
Consequently, as seen in (fl2l) . glj 2 ,} (0) = 0 at A 2 = 2, but 
not at A 2 = 2. 

The general result m also reveals that the classical 
Cauchy-Schwartz inequality for photon intensities is al¬ 
ways violated in the quantum case, i.e., 



(14) 
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Namely, introducing the parameter e = g bb \o)/g^a ( 0 ) 
the violation of the inequality requires [—(?ia' ) ( 0 )](l — 
y/e) 2 < 1/n which always applies since g^(0),n > 0 . 
Accordingly, emission of photons from the cavities occurs 
in a correlated way for all driving strengths and photon 
occupations. In the next Section we ascribe to the indi¬ 
vidual photon states in the cavities respective amplitudes 
(energies) and phases. One then sees that these states are 
correlated through their phases due to the simultaneous 
creation process in the transfer of a single Cooper pair. 


IV. SEMI-CLASSICAL DYNAMICS 


For large photon occupancy in the cavities a semiclas- 
sical type of approximation applies. The simplest semi- 
classical description of the dynamics of the system is ob¬ 
tained from the equations of motion for (a) and ( b ) which 
follow from Eq. 0, making the replacements (a) = a, 
( b ) = /? and treating expectation values of products of 
operators as products of expectation values. Hence we 
find 


a 




(*•' + £) 


a - 


iEj 

2A b Ej 


Ji(2A b \P\) x 


J 2 (2A 0 |a|). 


a 2 /3 


- J 0 (2A a |a|)( 


(15) 


1 


- ( i~5 {b ) + -)(3 + 


iEj 
2A a E c j 


Ji(2A a |a|) x 


J 2 ( 2 A b 


/ 3 2 a 


- Jo( 2 A & | 


a 


(16) 


where 8^’^ = ^( a > 6 ) / y/jajb- Obtained in this way, the 
factors of e ^ A “ +A ^/ 2 embodied in Ej that appear in 
these equations are accidental: they would not be present 
if we had instead chosen to use a symmetric ordering for 
the operators when deriving the Hamiltonian. However, 
Eqs. m and (1161) would also arise from a simple-minded 
ansatz in which we assumed that the density operator 
of the system is just a product of the coherent states 
p(t ) = \a(t))(a(t)\ ® |/3(/)}(/3(/)|, in this approximation 
the factors of e( A “ +A A / 2 would arise naturally. 

Using amplitude-phase coordinates for the two oscilla¬ 
tors, a = Ae~ l ^ a and /3 = Be _! ^, and introducing the 
total and relative phase variables = (f> a ±(j) bl Eqs. m 
and (fT 6 l) take the form 


A = 

B = 

t = 

r = 


E £ J 1 (2A b H)J 1 (2A a A) 

2 Ej 2A a A b A y ’ 

1 

~Yr B 


Ej J 1 (2A a A)J 1 (2A b H) 


E c , 


2A„A h B 


J ^L^. a L^. b l 

<5 (+) +F + (A, B) cos£ + 

8+ E_(A, B) cos£ + , 


sin(e + )(18) 

(19) 

( 20 ) 


where we used the Bessel function identity, 0 / 2 ( 2 :) + 
Jo (z) = 2Ji(z)/z, and have defined ± 5^. 


Further, 

E±(AB) = sSj ( ' Jl ( iy'’ [ J o(2A a A) - M2A a A)] 
± ' h{ \ a a B A) [Jo(2A b B) - J 2 {2A b B)] S j (21) 


with the property F + (—A,B) = F_(A,B) and 

F + (A, —B) = — F-(A, B). The behavior of the system is 
determined by the fixed points of the amplitudes A 0 ,B 0 
and the total phase . Since the relative phase does not 
appear on the righthand side of any of these equations 
its fixed point value is arbitrary. For simplicity, we con¬ 
centrate on the on-resonance case 8^ = 8^ = 0 in our 
analysis. 

The amplitude equations lead to the fixed point con¬ 
ditions A 0 = B 0 = 0 or 


rA a A b E c jA 2 

EjJi(2A b B 0 )Ji(2A a A 0 ) 

A a A b EjB 2 

rEjJi(2A b Bo)Ji(2A a Ao) 


( 22 ) 

(23) 


The second equality in Eq. [23] leads to the energy balance 
condition B 0 = rA 0 . From the equation for £ + , we see 
that fixed points arise when either cos £q“ = 0 or 


F + (A o ,B o )=0. (24) 


This latter condition is independent of Ej and hence 
leads to a locking of the amplitudes at particular val¬ 
ues as a function of Ej , something which is an important 
characteristic of the dynamics in the single-oscillator sys¬ 
tem—. For symmetric oscillators (r = 1 and A 0 = A b ) 
F + = 0 implies J[(z) = 0 with z = 2A a Ao = 2A b Bo 
which has a first solution at z = 1.841—. 

Thus we identify three possible fixed points for the sys¬ 
tem: a zero-amplitude one, one given by the conditions 
cos£ + = 0 and (from Eq. (1231) 1 


_ rAlA a A b E c j _ 

EjJi(2A b rAo)Ji(2A a Ao) 


= ± 1 , 


(25) 


and a third solution for which the amplitudes lock to 
values where Eq. (1M1) is satisfied (together with the con¬ 
dition B = rA) and the total phase is be given by Eq. 

m- 

We can look for small amplitude solutions to Eq. [25] 
(A b rA 0 , A a A Q < 1) by expanding the Bessel functions 
and retaining the lowest order terms in Aq , 


Aq = 



(26) 


Thus we see that a non-zero amplitude solution only ex¬ 
ists for Ej > Ej. Thus Ej has a simple physical inter¬ 
pretation: it is the value of Ej at which the oscillators 
reach the threshold for non-zero amplitude oscillations. 
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Taking into account the stability of the fixed points, we 
find that as Ej is increased from zero the amplitudes re¬ 
main zero until the system reaches threshold at Ej = Ej , 
after which the amplitudes grow smoothly according to 
Eq. (1251) with the global phase locked to £q~ = 7 t/2 . For 
a sufficiently large Ej , which we define as Ej 2 , a bifur¬ 
cation occurs as the amplitudes become large enough to 
satisfy Eq. (liMl) and the amplitudes then lock, becoming 
independent of Ej. 

In the next two sections we will examine the quantum 
dynamics of the system in the below and above threshold 
regimes. 


V. SUB-THRESHOLD DYNAMICS 

In the sub-threshold regime ( Ej < Ej) the semi- 
classical fixed points have zero amplitude (A = B = 0). 
In this case we can gain some insight into the behavior 
of the system by approximating the Hamiltonian of the 
system by its lowest order terms, i.e. setting Hrwa = 
h rwa t see Ec l- 0. an approach which is equivalent to 
analysing the linear fluctuations about the semi-classical 
fixed points. 

When this approximation is made the Hamiltonian is 
quadratic and the equations of motion for the moments 
take a rather simple form. Solving these equations, we 
find in the steady-state 



(a) = ( b ) = ( ab t) = 0 . (29) 


We note in passing that the result for n a reduces to the 
one derived in © in leading order in Ej/Ej. 

Simplified in this way, the linearized description leads 
to a Gaussian steady-state Wigner function which takes 
the for m 36 ' 37 


W a , b (a,i3) 


g-[(116+1/2) [a| 2 + (ra„ + l/2)|,3| 2 +/ia£+Ai*a*/3*]/C 


7T 2 C 


(30) 


where C = [(n a + 1/2)(rif, + 1/2) — \n\ 2 ] and p* = 
— (ab). This is a mixed state which combines two-mode 
squeezing and thermal-like fluctuations!!. The Wigner 
function of the individual oscillators is obtained by inte¬ 
grating over the phase space of the other one leading in 
either case to a thermal distribution. Thus for oscillator 
a, for example, we have 


W a (a) = 


n(n a + 1/2) 


exp 


(' n a + 1 / 2 ) _ 


(31) 




FIG. 3: (Color online) (a) Average occupation, ( n ) = n a = 
rib, (b) Fano factor, F = F a = Fb, as a function of Ej/Ej 
for symmetric oscillators. The full curves are the linearized 
results and the other curves are for A = 0.1 (dashed curves), 
A = 0.3 (dotted curves) and A = 0.6 (dash-dotted curves). 


The full behavior of the average energy of oscilla¬ 
tor a , n a , obtained by solving the master equation 
numerically 3 —, is shown in Fig. QT| for symmetric oscilla¬ 
tors (r = 1, A = A 0 = Ab). The divergence in n a which 
the linearized analysis predicts for Ej Ej [Eq. (|27|) [ 
never occurs in the full quantum problem as higher or¬ 
der terms in the RWA Hamiltonian always saturate the 
energy gain. As A is increased the saturation occurs at 
progressively lower values of the photon number whilst 
the range of Ej/Ej values for which the linearized cal¬ 
culation is accurate becomes smaller and smaller. 

The fluctuations in the energy of the oscillators, de¬ 
scribed by the Fano factors F a ^ J7J) change rather more 
dramatically with A a . The thermal Wigner function ob¬ 
tained from the linearized calculation [Eq. (13 1 fl ] predicts 
the simple relationship between Fano factor and photon 
number associated with thermal states, F a ^) = n a ^) + 1 , 
leading to growth in F a ^) as Fj/Ej increases and again 
there is a divergence at threshold. For small values of 
A, the full quantum dynamics follows a similar pattern 
though with saturation in F a ^) at the threshold leading 
to a peak rather than a divergence. In contrast, for larger 
A values the behavior is completely different: the value 
of F a (p) drops monotonically as Ej/Ej is increased and 
its behavior contains no signature of the threshold at Ej. 

The change in the behavior of F a as A a is increased is 
reminiscent of quantum optical systems like the laser—, 
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which in the ‘thermodynamic’ limit of weak atom-photon 
couplings display clear thresholds (accompanied by a sig¬ 
nature peak in the Fano factor) whose properties can be 
understood in terms of an analogy with classical phase 
transitions, but which for sufficiently strong couplings be¬ 
have quite differently without clear signatures of a thresh¬ 
old—^. 


VI. DYNAMICS ABOVE THRESHOLD 

Above threshold the oscillators become strongly ex¬ 
cited though this does not mean that their states become 
classical. As in the case of the single-oscillator system—, 
strong number squeezing (marked by a Fano factor below 
unity) occurs even at large average occupation numbers. 
As in the sub-threshold regime, the behavior of the sys¬ 
tem in the limit of very small zero-point fluctuations, 
A q , Ab -C 1, can be captured within an approximate de¬ 
scription which linearizes about the semi-classical fixed 
points of the system, but for larger zero-point fluctua¬ 
tions numerical solution of the quantum master equa¬ 
tion becomes essential. We start by exploring the general 
properties of the steady-states of the individual oscilla¬ 
tors in the above-threshold regime for symmetric oscil¬ 
lators and the role played by the size of the zero-point 
fluctuations before going on to examine how asymmetry 
alters the behavior. 


A. Symmetric Oscillators 

For symmetric oscillators (r = 1, A Q = A{, = A) 
the steady-state properties of the two oscillators must be 
the same and there is a very simple scaling to the semi- 
classical fixed point amplitudes obtained in Sec. IIVI the 
value of 2A a A 0 is a function of just Ej/Ej , see (E3l) . 
This scaling provides a convenient way of comparing the 
average oscillator occupation n = n a = rib (obtained by 
solving the master equation numerically) for different val¬ 
ues of A with the semi-classical prediction, as shown in 
Fig. UJi. We solved the master equation using standard 
numerical methods—; for smaller values of A we carried 
out quantum trajectory simulations, whilst for larger A 
we were able to solve for the steady-state of the master 
equation directly because the state-space required was 
rather smaller. Indeed, the strong suppression in the 
magnitude of the oscillator occupation number as A is 
increased (there is a reduction by a factor ~ 100 in going 
from A = 0.1 to A = 0.6) is the most significant feature 
in Fig. [4^, which is captured by the 4A 2 scaling. 

Figure[3Ji also shows that the semi-classical amplitudes 
provides a very good description of the oscillator occu¬ 
pations for A < 1. For A = 0.1 we see that there are 
small deviations from the semi-classical predictions which 
become apparent just above threshold and near the bi¬ 
furcation that occurs at E c2 = 2.5 Ej. As the size of the 
zero-point fluctuations is increased, these small devia- 



FIG. 4: (Color online) Comparison of the oscillator occupa¬ 
tion numbers (a) and Fano factor (b) obtained from numeri¬ 
cal solution of the quantum master equation for A = 0.1, 0.3 
and 0.6 with corresponding semi-classical calculations over 
the range Ej < Ej < E 1 } 2 = 2.5 Ej. In (a) both the semi- 
classical oscillator energy, Aq, and occupation number, n a are 
scaled by 4A 2 . 


tions grow much larger and spread out over a much wider 
range of Ej/Ej values. Nevertheless, the semi-classical 
amplitude continues to provide a useful estimate of the 
full quantum results even for A = 0.6. 

We now turn to the fluctuations in the occupation 
numbers of the oscillators, described by the single mode 
Fano factors, F = F a = Ft. The value of F decreases 
progressively the further above threshold we go as shown 
in Fig. [4jr. For very small A, F is strongly elevated 
close to threshold (the other side of the peak in F seen 
below threshold), but decreases rapidly with increasing 
Ej/Ej leading to substantial number state squeezing 
with F ~ 0.5 before the bifurcation at E c2 . For larger 
A values there is no peak around threshold and F < 1 
throughout though the lowest values are slightly larger 
than those obtained for very small A. 

The simple semi-classical analysis in Sec. IIVI can be 
extended to describe fluctuations (at least to lowest or¬ 
der) by essentially adding a noise term to the equations 
of motion for the amplitudes, Eqs. m and csi) , so 
that they become Langevin equations. Formally, such 
Langevin equations can be derived within the frame¬ 
work of an approximate semi-classical approach known 
as the truncated Wigner approximation, as we show in 




Appendix [A] We again make the change to amplitude- 
phase variables and then linearize about the fixed point 
values to obtain expressions for the amplitude fluctua¬ 
tions ( SA 2 ) = ((A — Ao) 2 ) which can be related to the 
Fano factor in a simple way F a ~ 4((5A 2 ) (details of the 
calculation are provided in Appendix [All. 

The comparison of the semi-classical and quantum cal¬ 
culations of the Fano factor shown in Fig. |3b shows that 
the semi-classical Fano factor, which is a function of 
Ej/Ej alone in the symmetric case, can be thought of 
as giving the low-A limit. As A is increased the devi¬ 
ations from the semi-classical value get stronger around 
threshold and the bifurcation at E c2 = 2.5Ej as well as 
spreading over a wider range of Ej / Ej in much the same 
way as for the oscillator occupation. Note that the semi- 
classical calculation predicts a Fano factor which tends 
to 0.5 as the system tends to the bifurcation, Ej —> E c2 . 
This matches the lowest Fano factors found for the one- 
oscillator system which occurs as the system tends to¬ 
wards an above-threshold bifurcation at the 2-photon res- 

2D 

onance—. 


B. Asymmetric oscillators 

We now consider what happens when the oscillators 
are no longer entirely symmetric. We start by considering 
the case where the zero-point fluctuations of the modes 
remain the same (A = A a = Af ,), but the damping rates 
are different r yf 1 and then go on to consider the general 
case where A 0 ^ Af, and r yf 1. 

The effect of asymmetric damping on the average oc¬ 
cupation numbers of the oscillator (shown in Fig. [5]), is 
twofold with both effects following from the underlying 
semi-classical dynamics discussed in Sec. IIVI Firstly, the 
bifurcation which occurs at E c2 is pushed to larger values 
of Ej. Secondly, the average energies of the modes be¬ 
come unequal in proportion to the underlying asymmetry 
in the damping, rib = r 2 n 0 . 

Figure |6] shows the effect of asymmetric damping on 
the occupation number fluctuations for different values of 
A. What is striking here is that the fluctuations become 
asymmetric and the Fano factor becomes significantly 
lower than 0.5 in the less damped oscillator. The low¬ 
est values of F are achieved well-above threshold, close 
to the bifurcation at Ej 2 for small-A, though for larger 
A values the minimum A is at a lower value of Ej as 
the increase in F associated with the bifurcation starts 
to occur at progressively smaller values of Ej/Ej as A 
is increased. Above the bifurcation the value of F settles 
down to a steady, but rather higher value. 

The semi-classical calculation predicts a minimum 
value of F ~ 0.1 for the small-A limit when r = 1/3, sub¬ 
stantially lower than any of the Fano factors predicted for 
the single-oscillator system—, and this value continues to 
decrease for smaller r. This suggests that the asymmetric 
two-oscillator system may provide a very effective route 
to preparing a particular mode in a strongly non-classical 



0 2 4 6 8 10 12 14 


FIG. 5: (Color online) Steady-state occupations n a and rib 
(full lines) compared with the classical values of Aq and Bq at 
the stable fixed points (dashed lines) for (a) A = 0.3 (b) A = 
0.6. Results are shown for r = 1, 1/2 and r = 1/3 in each case. 
Note that the semi-classical amplitudes are zero for Ej < Ej. 
The above-threshold bifurcation occurs at Ej 2 /Ej = 2.5, 4.0 
and 8.7 for r = 1, 1/2 and r = 1/3, respectively. 


state at large photon numbers. As F —> 0 the state of 
the oscillator must eventually become a pure Fock state 
and so one naturally expects to find negative features in 
the Wigner function for very small values of F. However, 
the presence of negative regions in a Wigner function is 
not simply a function of F , but also the average occupa¬ 
tion number (n): as one goes to larger average oscillator 
occupation numbers, smaller and smaller values of F are 
required to form negative regions. Figure [7] illustrates 
this by showing examples of the Wigner functions for 
A = 0.3 and A = 0.6 with r = 1/3 and Ej/Ej = 6 where 
F ~ 0.2 in both cases (see Fig. [6]). For A = 0.6 there 
is strong evidence of negativity in the Wigner function 
whilst it is almost washed out for A = 0.3 since although 
the Fano factors are very similar, the latter has a much 
higher average occupation number. 

Finally, we examine the behavior in the regime where 
A a Af,. Figure [8] shows examples of the behavior of the 
occupation numbers and Fano factors of the two oscilla¬ 
tors in this case. Interestingly for r = 1 whilst energy 
balance means that n a = rib, the fluctuations in the two 
modes are no longer the same. When r / 1 the occu¬ 
pation numbers of the two oscillators spilt according to 
the usual relation, rib = r 2 n a and the fluctuations be¬ 
come even more asymmetric. Indeed, the minimum val¬ 
ues of the Fano factors, are lower than those in the cor- 
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FIG. 6: (Color online) Steady-state Fano factors of the modes 
(a) calculated semi-classically (small-A limit) and calculated 
numerically using the master equation for (b) A = 0.3 and 
(c) A = 0.6. In each case results are shown for r = 1 (full 
lines), r = 1/2 (dashed lines) and r = 1/3 dotted lines. For 
r = 1/2 and r = 1/3 the upper curves are for oscillator b 
and the lower ones for oscillator a. The bifurcation occurs 
at Ej 2 /Ej = 2.5, 4.0 and 8.7 for r = 1, 1/2 and r = 1/3, 
respectively. Note that the semi-classical results in (a) are for 
Ej < Ej < Ej 2 whilst (b) and (c) cover a broader range of 
Ej values. 


responding cases where A a + A& takes the same value, 
but A a = A b . 


VII. DISCUSSION AND CONCLUSIONS 

We have analyzed the quantum dynamics of two 
electromagnetic oscillators coupled to a voltage biased 
Josephson junction. We considered the case where the 


FIG. 7: (Color online) Wigner function of oscillator a for 
r = 1/3, Ej/Ej = 6 and (a) A = 0.3 (b) A = 0.6. Negative 
regions are apparent in both cases, though more strongly in 
(b). The Fano factors associated with the states are F a = 0.19 
(a) and F a = 0.22 (b). 


voltage across the junction was tuned so that the energy 
lost by a Cooper pair crossing the circuit matches the 
sum of the photon energies of the two oscillators. In this 
regime the oscillators are pumped by the flow of Cooper 
pairs and can become strongly excited. Using a rotat¬ 
ing wave approximation, we derived an effective time- 
independent Hamiltonian for the system and explored the 
behavior it gives rise to under a wide range of conditions 
using a mixture of numerical and analytic approaches to 
solve the master equation. We use a perturbative ap¬ 
proach to obtain analytic results for the regime where 
the occupation of the oscillators is low while in the oppo¬ 
site regime of large occupation numbers a semi-classical 
approach provides an effective description. 

The steady states of the oscillators display signatures 
of non-classical behavior over a very wide range of condi¬ 
tions with sub-Poissonian photon statistics found in both 
the low and high occupancy regimes. The strength of the 
zero-point fluctuations in the oscillators, A a ( b ), plays an 
important role: as these are increased the overall excita¬ 
tion level of the oscillators tends to move towards lower 
photon numbers whilst the signatures of non-classicality 
are enhanced. The ratio of the damping rates of the two 
cavities, described by r = \Jr a /r bl also has an interesting 
effect on the behavior of the system. The photon num¬ 
bers in the two oscillators are related in a simple way, 
n b = r 2 n a , as one would expect. However, the quantum 
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FIG. 8: (Color online) Steady-state occupation numbers, n, 
(a) and Fano factors, F, (b) of the oscillators for A a = 0.4, 
A ( , — 0.2 with r = 1 and 1/2. In (a) the semi-classical pre¬ 
dictions are shown as a dotted line and the numerical results 
as full line in each case. 


fluctuations (e.g. measured by the Fano factors F a ( b )) also 
become unequal in the asymmetric case, r / 1. Indeed 
we find that the Fano factor in the less-damped oscilla¬ 
tor can become low enough to lead to significant negative 
regions in the corresponding Wigner function. 


Strong correlations between the two oscillators are to 
be expected in the regime we consider given the fact that 
the tunnelling Cooper-pairs excite photons in each of 
the two oscillators simultaneously. The violation of the 
classical Cauchy-Schwarz inequality for the photons in 
the two oscillators, g ab , indicates that the corresponding 
two-mode states are non-classical. It would be natural 
to also investigate the entanglement between the two os¬ 
cillators. However, this is complicated by the fact that 
in practice local voltage fluctuations, even when weak, 
would be expected to have a very strong influence on 
phase dependent correlation functions such as ( ab ) which 
can be important in determining the level of entangle¬ 
ment. This is in contrast to the observables such as pho¬ 
ton occupation numbers and correlation functions which 
we have focussed on here which, as remarked in Sec. Ill Bl 
are expected to be only very weakly affected. We plan to 
address the issue of inter-oscillator entanglement in a fu¬ 
ture work using a form of the master equation where the 
effects of voltage fluctuations are explicitly included^!. 
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Appendix A: Semi-classical calculation of 
above-threshold fluctuations 


We can gain useful insights into the dynamics by 
extending our semi-classical analysis to include quan¬ 
tum fluctuations using a truncated Wigner approxima¬ 
tion (TWA) 31 i 40 . The TWA leads to an approximate 
equation of motion for the Wigner function of the sys¬ 
tem, W (a, /3), in which third-order and higher derivatives 
are neglected. Dropping higher-order derivatives leads 
to a Fokker-Planck equation from which we obtain^! 
Langevin equations for the phase space variables a, j3 
of the form (for the on-resonance case) 


a 


$ 


iEj 


r 

- ol , . „ 

2 2A h E' 


■Ji(2A b 


b&j 

J 2 (2A a H)^-- J 0 (2A a |a|)^ 


|a| 

—L p + lEj 
2 r 2A a Ej 


Ji(2A 0 |a|) x 


J 2 ( 2A fc 


/3 2 a 


a 


a 


- Jo(2A b |/3|) — 


a 


(Al) 

+ Va(t) 

(A2) 

+ vp(t). 


The noise terms r] a (/ 3 )(t) have zero means and the only 
non-zero second moments are given by 

= r -8{t-t') (A3) 

M^*(*')> = (A4) 

Apart from the noise terms, the equations of motion take 
the same forint as those derived in Sec. lIVI [Eas. (fl5l) and 

mi 

We proceed by changing to amplitude and phase vari¬ 
ables and then linearizing about the fixed point values, 
i.e. working to first order in 6A = A — Ao, SB = B — Bq 
and c>£ + = £ + — with A 0 , B 0: ^ the fixed point 
values. For the fixed point just above threshold the am¬ 
plitude and phase fluctuations become decoupled and on- 
resonance we And 


(5A 

SB 


-r Q V,b) \ f \ , ( \ 

h{b,a) — r b ) \ SB J \Vb ) 


(A5) 
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where 


r 


a 


r b 


a,b) 


- + ( ^ aEj ^ X 

2 + V2A b E c j) 

Ji(2A b B Q ) [J\(2A a A 0 ) + ^3(2A a ^4 0 )] 

J_ / A^\ 

2r^\2A a E c j) 

Jr(2A a A 0 ) [Ji(2A b B 0 ) + J 3 (2A b B 0 )} 


(A6) 

(A7) 


(^r) [Jo(2A b S 0 ) - J 2 (2A b B 
[J 0 (2A a A 0 ) +J 2 (2A a A 0 )] 


0 


)] X (A8) 


and a corresponding expression for h^ ba y The noise 
terms obey the correlation functions 


{VA{t)r] A {t')) = 7 8 9 10 11 12 --6(t-t’) (A9) 

(■ riB{t)vB(t ')} = l~5{t-t'). (A10) 


Using Eq. (1A5I) we obtain the steady-state variances 


(<L4 2 ) 
(SB 2 ) 
{SASB) 


+ ^!l(6A5B) (All) 

-L- + ^A±(6A5B) (A12) 

8rT b T b 

h ia , b )T a /r + ft (b , n) r b r 

8(r Q + r b )(r a r b - /i( 0 ,6)fy&,a)) 


Recalling that a and /? are phase space variables of 
a Wigner function, we can connect these variances to 
quantum averages: (A 1 2 ) = (a)a) + 1/2 and (A 4 ) = 
((a^a) 2 ) + {ata) + 1/2. For fixed points where A 0 3> 1, 
corrections of order Aq 2 can be neglected, leading to the 
simple result, 


(A 4 ) - (A 2 ) 2 - 1/4 

(A 2 ) - 1/2 

4Ag(5A 2 ) + 2{8A 2 ) 2 — 1/4 
~ An + {SA 2 ) - 1/2 
=* 4(<5A 2 ), 


(A14) 

(A15) 

(A16) 


and there is of course a corresponding relation for F b . 
The Langevin equation for <5£ + takes the form 

= -F + {Ao,B 0 )Sf + Vi+, (A17) 

where (r]£+ {t)rj^+ (t')) = 2D6(t — t') with 2D = r/(4Ag) + 
l/(4r\Bg). Hence we find 

m+) 2 )=D/F + (A 0 ,B 0 ). (A18) 


Note that as the system approaches the bifurcation at 
Ej = Ej 2 , F + (A 0 ,B 0 ) —> 0 implying that the total phase 
fluctuations within this linearized approach diverge. 
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